%matplotlib notebook
import numpy as np
import scipy.signal
import scipy.fft as sfft
import matplotlib.pylab as plt
from matplotlib import animation

from IPython.display import YouTubeVideo, HTML, Audio
from bokeh.layouts import column, row
from bokeh.models import CustomJS, ColumnDataSource, Slider
from bokeh.plotting import Figure, show, output_notebook
output_notebook()
Loading BokehJS ...

2. Sistemas para el procesamiento de señales

2.1. Definición de sistema

Hasta ahora hemos realizado análisis de señales, es decir el estudio de las señales y sus propiedades en el dominio del tiempo y frecuencia

En esta unidad nos enfocaremos en el procesamiento de señales, es decir el diseño de sistemas que procesan una señal de entrada y producen una señal de salida

../../_images/system.png

Usaremos

  • \(x[n]\) para denotar la señal (discreta) de entrada y \(X[k]\) su espectro

  • \(y[n]\) para denotar la señal (discreta) de salida e \(Y[k]\) su espectro

2.1.1. Ejemplos de sistemas

Utilizando sistemas podemos modificar una señal para mejorar su calidad o remover efectos indeseados

  • Un sistema para reducir el ruido de una señal de electroencefalograma (EEG)

../../_images/system-denoise-eeg.png
  • Un sistema para mejorar una imagen fuera de foco (sharpening)

../../_images/system-sharpen.jpg
  • Un sistema para eliminar el eco de un audio

../../_images/system-echo.png

En lo que sigue realizaremos una clasificación de los sistemas según sus propiedades

2.2. Sistemas sin memoria

Diremos que un sistema \(\Phi\) es un sistema sin memoria si

\[ y[n] = \Phi(x[n]), \]

es decir que la salida del sistema en un instante \(n\) dado depende solo de la entrada en ese mismo instante

Veamos algunos ejemplos

2.2.1. Sistema atenuador/amplificador ideal

\[ y[n] = A x[n], \]

donde \(A>0\) se llama ganancia

Este sistema puede atenuar la entrada si \(0<A<1\) y amplificar si \(A>1\)

2.2.2. Sistema saturador (clamp)

\[\begin{split} y[n] = \begin{cases} B &x[n] > B \\x[n] & x[n] \in [-B, B]\\ -B & x[n] < -B\end{cases} \end{split}\]

Este sistema limita los valores de la señal de entrada en un rango fijo

2.2.3. Sistema rectificador

\[ y[n] = | x[n] | \]

Este sistema eliminar la parte negativa de la señal de entrada

2.3. Sistema Lineal

Diremos que un sistema \(\Phi\) es lineal si cumple con las siguientes propiedades

2.3.1. Homogeneidad

Un cambio en la amplitud de la entrada produce un cambio equivalente en la salida

\[ \Phi(cx[n]) = c \Phi(x[n]) = c y[n] \]

2.3.2. Aditividad

Señales que se suman en la entrada producen señales que se suman en la salida

\[ \Phi(x_1[n] + x_2[n]) = \Phi(x_1[n]) + \Phi(x_2[n]) = y_1[n] + y_2[n] \]

Es decir que las señales pasan por el sistema sin interactuar entre ellas

2.3.3. Otras propiedades de los sistemas lineales

Producto de las propiedades anteriores se tiene que una cascada de sistemas lineales forma un sistema lineal equivalente

Y la cascada de sistemas es conmutativa: El orden de los sistemas en la cascada no altera el resultado final

../../_images/system-conmu.png

Los sistemas lineales también cumplen el Principio de superposición

  1. Si descomponemos una señal en \(M\) componentes: \(x[n] = x_1[n] + x_2[n] + \ldots + x_M[n]\)

  2. Y aplicamos un sistema lineal a cada componente: \(y_j[n] = \Phi(x_j[n])\)

  3. Podemos recuperar la salida total usando aditividad: \(y_1[n] + y_2[n] + \ldots + y_M[n] = y[n]\)

../../_images/system-superpos.png

2.4. Sistemas con memoria

Un sistema \(\Phi\) es un sistema con memoria si su salida actual depende sólo de la entrada actual, las entradas anteriores o las salidas anteriores

\[\begin{split} \begin{align} y[n] = \Phi(x[n], & x[n-1], x[n-2], \ldots, x[0], \\ \nonumber & y[n-1], y[n-2], \ldots, y[0]) \nonumber \end{align} \end{split}\]

esto también se conoce como sistema causal

Un sistema con memoria no-causal usa entradas futuras (es decir \(x[n+1]\), \(x[n+2]\), …) y por ende solo se puede implementar de forma offline, es decir una vez que sea ha observado toda la señal

Veamos algunos ejemplos de sistemas con memoria

2.4.1. Sistema con un retardo (delay)

Definido como

\[ y[n] = x[n-m], \]

donde

  • la salida depende solo de «una» entrada anterior

  • el valor de m define que tan «antigua» es la entrada pasada

El delay o retarno no afecta la amplitud de los componentes frecuenciales de la señal pero si su fase, como muestra la siguiente animación

%%capture

fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m)) 
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))

def update(m):
    ax[0].cla(); ax[0].plot(n, x(m));
    X = sfft.fftshift(sfft.fft(x(m)))
    Xm = np.absolute(X); Xp = np.angle(X)
    # Espectro de magnitud:
    ax[1].cla(); ax[1].plot(f, Xm); 
    # Espectro de fase enmascarado con el espectro de magnitud
    ax[2].cla(); ax[2].plot(f, Xm*Xp/np.amax(Xm)); 
    ax[2].set_ylim([-np.pi, np.pi])
    angle_delay = Xp[np.argmax(Xm)]
    ax[2].set_title("%0.4f [rad], %0.0f [deg]" % (angle_delay, 180*angle_delay/np.pi))
    return ()

anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=True)
HTML(anim.to_html5_video())

2.4.2. Sistema reverberador o eco

Definido como

\[ y[n] = x[n] + A x[n-m], \]

donde

  • la salida depende de una entrada «pasada» y la entrada actual

  • la ganancia controla si el «eco» se atenua o amplifica

Al contrario del sistema anterior, el eco si puede modificar el espectro de amplitud.

Notemos el efecto de interferencia constructiva y destructiva al modificar el retardo, como muestra la siguiente animación

%%capture
fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m)) 
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))
A = 1.
def update(m):
    y = x(0) + A*x(m)
    X = sfft.fftshift(sfft.fft(y))
    for ax_ in ax:
        ax_.cla()
    ax[0].plot(n, x(0), label='x[n]')
    ax[0].plot(n, A*x(m), label='A*x[n-m]')
    ax[0].legend()
    ax[1].plot(n, y); 
    ax[1].set_ylim([-A-1.1, A+1.1])
    ax[2].plot(f, np.absolute(X)); 
    ax[2].set_ylim([-20, (1+A)*len(n)/2 + 20])

anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=False)
HTML(anim.to_html5_video())

Ejemplo de interferencia destructiva en una onda mecánica: https://www.youtube.com/watch?v=IU8xeJlJ0mk

2.4.3. Sistemas con múltiples ecos

Pueden combinarse más retardos para hacer un sistema reverberante más complejo

Por ejemplo

\[ y[n] = x[n] + A_1 x[n-m_1] + A_2 x[n-m_2] + A_3 x[n-m_3] + \ldots, \]

Que se escucha como

Fs = 44100; 
n = np.arange(0, 4, step=1.0/Fs) 
x = lambda m: np.sin(2.0*np.pi*880*(n-m))*np.exp(-(n-m)**2/0.5**2)*np.heaviside(n-m, 0)
y = x(0) + 0.5*x(1.) + 0.25*x(2.) + 0.125*x(3.)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(n, y,  line_width=3)
p1.xaxis[0].axis_label = 'Tiempo [s]'
show(p1)
Audio(y, rate=Fs, normalize=False)

2.5. Sistema FIR

Generalizando el ejemplo de sistema lineal reverberante a \(L\) retardos llegamos a

\[\begin{split} \begin{align} y[n] &= h[0] x[n] + h[1] x[n-1] + h[2] x[n-2] + \ldots + h[L] x[n-L] \nonumber \\ &= \sum_{j=0}^{L} h[j] x[n-j] \nonumber \\ &= (h* x)[n] \nonumber \end{align} \end{split}\]

que se puede modelar como una convolución discreta entre \(h\) y \(x\)

Este sistema se conoce como

  • Sistema FIR (finite impulse response)

  • Sistema MA (moving average)

  • Sistema todo-zeros

y es de orden L (posee L+1 coeficientes)

2.5.1. Intepretación como media movil (MA)

El sistema FIR es equivalente a una media movil ponderada que se aplica sobre la entrada

Los coeficientes del sistema son los ponderadores

Por ejemplo sea un sistema de 3 coeficientes \(h[0]=a\), \(h[1]=b\) y \(h[2]=c\)

\[\begin{split} \begin{align} y[n] = (h*x)[n] &= \sum_{j=0}^{2} h[j] x[n-j] \nonumber \\ &= a x[n] + b x[n-1] + c x[n-2] \nonumber \end{align} \end{split}\]

donde cada salida se calcula a partir de

\[ \overbrace{x[0], x[1], x[2]}^{y[2]} , x[3], x[4], \ldots \]
\[ x[0], \overbrace{x[1], x[2] , x[3]}^{y[3]}, x[4], \ldots \]
\[ x[0], x[1], \overbrace{x[2] , x[3], x[4]}^{y[4]}, \ldots \]

Para obtener el valor de \(y[0]\) e \(y[1]\) se deben establecer «condiciones de borde», como por ejemplo \(x[-2] = x[-1]= 0\)

2.5.2. Ejemplo: Eliminando ruido blanco aditivo

Si tenemos una señal contaminada con ruido blanco aditivo podemos usar un sistema FIR promediador para «suavizar la contaminación»

La animación muestra un sistema FIR con \(L\) coeficientes idénticos e iguales a \(1/L\) que se convoluciona con la señal contaminada intentando recuperar la «señal limpia»

%%capture
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C) 
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
# Sistema
L = 10; h = np.ones(shape=(L,))/L; 
# Acumulator
y = np.zeros_like(n, dtype=np.float)

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
ax[0].plot(n, x, 'k.'); 
ylims = ax[0].get_ylim()

def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()

def update(m):
    ax[1].cla(); ax[2].cla()
    c = np.zeros_like(n, dtype=np.float); 
    c[m:m+L] = h
    ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+L])
    ax[2].plot(n, y, label='señal recuperada');  
    ax[2].plot(n, x_clean, 'g', label='señal limpia'); 
    ax[2].legend()
    ax[2].set_ylim(ylims)
    ax[2].plot([m, m], [ylims[0], ylims[1]], 'r--', alpha=0.5)
    return ()

anim = animation.FuncAnimation(fig, update,init_func=init, 
                               frames=100-len(h), interval=150, blit=True)
HTML(anim.to_html5_video())

2.5.3. Ejemplo: Encontrando cambios en una señal

Podemos usar un sistema «diferenciador» para detectar cuando una señal cambia de valor como se muestra en la siguiente animación

%%capture

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
n = np.arange(0, 100, step=1)
x = np.zeros_like(n, dtype=np.float)
x[20:] += 1.; x[40:] += 1.; x[80:] -= 1.;
ax[0].plot(n, x)
# System:
h = np.array([-0.5, 0.5])
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()

def update(m):
    c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
    ax[1].cla(); ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+len(h)])
    ax[2].cla(); ax[2].plot(n, y);  
    ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
    return ()

anim = animation.FuncAnimation(fig, update, init_func=init, 
                               frames=100-len(h), interval=150, blit=True)
HTML(anim.to_html5_video())

2.5.4. Ejemplo: Removiendo una tendencia

En el ejemplo siguiente tenemos una señal sinusoidal que nos interesa separar de una señal más lenta/suave denominada tendencia

Podemos usar un sistema FIR para eliminar la tendencia

%%capture

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
np.random.seed(0); 
n = np.arange(0, 100, step=1)
C = np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/30**2)
x = np.sin(2.0*np.pi*0.1*n) + 5*np.random.multivariate_normal(np.zeros_like(n), C)
ax[0].plot(n, x); ylims = ax[0].get_ylim()
# System:
L = 5; h = -np.ones(shape=(L,))/L; h[L//2] += 1
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()
    
def update(m):
    c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
    ax[1].cla(); ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+len(h)])
    ax[2].cla(); ax[2].plot(n, y);  
    ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
    return ()
    
anim = animation.FuncAnimation(fig, update, init_func=init, 
                               frames=100-len(h), interval=200, blit=True)
HTML(anim.to_html5_video())

2.6. Convolución con scipy

En la práctica, si queremos convolucionar una señal con nuestro sistema, podemos usar la función de scipy

scipy.signal.convolve(in1, # Señal de entrada
                      in2, # Coeficientes del sistema
                      mode='full', 
                      method='auto'  
                     )

donde el argumento method puede ser

  • direct: Realiza la convolución en el dominio del tiempo

  • fft: Realiza la convolución multiplicando los espectros

  • auto: Se decide automaticamente en base al largo de las señales

y el argumento mode indica donde se hace la convolución

Sea una señal \(x=[a,b,c]\) y un sistema \(h=[d, e]\)

  • Si uso mode=valid el resultado será \(y=[ad+be, bd+ce]\)

  • Si uso mode=same el resultado será \(y=[ae, ad+be, bd+ce]\), es decir se agregan ceros al principio de \(x\) tal que \(y\) sea del mismo largo que \(x\)

  • Si uso mode=full el resultado será \(y=[ae, ad+be, bd+ce, cd]\), es decir se agregan ceros al principio y al final de \(x\)

2.6.1. Ejemplo: Eliminando ruido v 2.0

Esta vez usamos scipy.signal.convolve y probamos dos sistemas FIR

  • coeficientes idénticos e iguales a \(1/L\) (rectangular)

  • coeficientes que decaen suavemente a cero (hamming)

from bokeh.palettes import Dark2_5 as palette
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C) 
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
L = 20; 

p = [Figure(plot_width=600, plot_height=230, toolbar_location="below") for k in range(3)]

p[0].scatter(n, x, color='black')
p[2].line(n, x_clean, line_width=2, legend_label = 'señal limpia'); 
for name, h, color in zip(['rect', 'hamming'],
                          [np.ones(shape=(L,)), scipy.signal.hamming(L)],
                          palette):

    h = h/np.sum(h)
    p[1].line(np.arange(L), h, color=color, legend_label=name, line_width=2)
    p[2].line(n, scipy.signal.convolve(x, h, mode='same', method='auto'), 
              color=color, legend_label=name, line_width=2); 
p[2].legend.location = "top_left"
show(column(p))